function testspec(sine_amp, sine_freq, noise_amp, window_fcn);

% TESTSPEC - test Lomb & FFT PSD methods using simulated data
%
%    TESTSPEC(AMP, FREQ, NOISE, WIN) calculates the PSD of sinusoids having
%    specified amplitudes and frequencies using Lomb's method and FFTs.  AMP 
%    is a vector of sinusoidal amplitudes, FREQ is the corresponding vector
%    of sinusoidal frequencies, and NOISE is a scalar to provide the standard
%    deviation of the Gaussian noise.  An optional string WIN provides the
%    name of a MATLAB windowing function ('hanning', etc.) to ensure that
%    windowing does not affect the results.
%
%    To simulate the varying length of data, a time base will created from a
%    random number of points and random duration.  Note that the resulting
%    sampling rate will not fall below 2*MAX(FREQ) to prevent aliasing. 

% By:   S.C. Molitor (smolitor@eng.utoledo.edu)
% Date: February 1, 2001

% validate arguments

if ((nargin < 3) | (nargin > 4))
   msgbox('Invalid number of arguments', 'TESTSPEC Error', 'warn');
   return
elseif (~isnumeric(sine_amp) | isempty(sine_amp))
   msgbox('AMP must be a numeric vector', 'TESTSPEC Error', 'warn');
   return
elseif (~isnumeric(sine_freq) | isempty(sine_freq))
   msgbox('FREQ must be a numeric vector', 'TESTSPEC Error', 'warn');
   return
elseif (prod(size(sine_amp)) ~= prod(size(sine_freq)))
   msgbox('AMP and FREQ must have the same number of elements', 'TESTSPEC Error', 'warn');
   return
elseif (~isnumeric(noise_amp) | (prod(size(noise_amp)) ~= 1))
   msgbox('NOISE must be a scalar', 'TESTSPEC Error', 'warn');
   return
elseif (nargin == 3)
   window_fcn = '';
elseif (~ischar(window_fcn))
   msgbox('WIN must be a string', 'TESTSPEC Error', 'warn');
   return
end

% generate time base from random number of points & sampling rate

num_pts = round(200 + 800 * rand(1, 1));		% NPTS between 200 - 1000
samp_freq = max(sine_freq) * (20 + 80 * rand(1, 1));	% SAMP FREQ between 20 - 100 * Fmax
x_time = [1 : num_pts] / samp_freq;

% generate data with noise & sinusoids
% sinusoids have specified amp/freq, random phases

y_data = noise_amp * randn(size(x_time));
for i = 1 : length(sine_amp)
   y_data = y_data + sine_amp(i) * sin(2 * pi * (sine_freq(i) * x_time - rand(1, 1)));
end
subplot('Position', [0.1, 0.55, 0.8, 0.4]);
plot(x_time, y_data, 'b-');

% generate frequencies over which to evaluate PSDs
% highest frequency is 10*MAX(FREQ)

x_freq = 10 * max(sine_freq) * [1 : num_pts] / num_pts;

% window data if specified

if (isempty(window_fcn))
   win_data = ones(size(x_time));
else
   win_data = feval(window_fcn, length(x_time));
   win_data = win_data';
end
win_mag = sum(win_data) / length(win_data);
y_data = y_data .* win_data;
hold on
plot(x_time, y_data, 'r-');
hold off

% evaluate PSDs using FFT

n_fft = floor(length(y_data) / 2);
y_fft = fft(y_data);
x_fft = (samp_freq / 2) * [0 : n_fft] / n_fft;
p_fft = sqrt(abs(y_fft(1 : n_fft + 1)) .^ 2) / (n_fft * win_mag);
i_fft = find(x_fft <= max(x_freq));

% evaluate PSDs using Lomb

p_lomb = lombpsd_scm(x_time, y_data, x_freq);
p_lomb = sqrt(p_lomb / n_fft) / win_mag;
subplot('Position', [0.1, 0.1, 0.8, 0.4]);
plot(x_freq, p_lomb, 'r-');
hold on
plot(x_fft(i_fft), p_fft(i_fft), 'bx');

% re-evaluate Lomb PSD at FFT frequencies

p_lomb2 = lombpsd_scm(x_time, y_data, x_fft(i_fft));
p_lomb2 = sqrt(p_lomb2 / n_fft) / win_mag;
plot(x_fft(i_fft), p_lomb2, 'go');
hold off
return
